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Abstract. The dynamic model for large-eddy simulation (LES) of turbulent flows requires test fil- 
tering the resolved velocity fields in order to determine model coefficients. Ifowever, test filtering is 
costly to perform in large-eddy simulation of complex geometry flows, especially on unstructured 
grids. The objective of this work is to develop and test an approximate but less costly dynamic 
procedure which does not require test filtering. The proposed method is based on Taylor series 
expansions of the resolved velocity fields. Accuracy is governed by the derivative schemes used 
in the calculation and the number of terms considered in the approximation to the test filtering 
operator. The expansion is developed up to fourth order, and results are tested a priori based 
on direct numerical simulation data of forced isotropic turbulence in the context of the dynamic 
Smagorinsky model. The tests compare the dynamic Smagorinsky coefficient obtained from fil- 
tering with those obtained from application of the Taylor series expansion. They show that the 
expansion up to second order provides a reasonable approximation to the true dynamic coefficient 
(with errors on the order of about 5% for Cj), but that including higher-order terms does not nec- 
essarily lead to improvements in the results due to inherent limitations in accurately evaluating 
high-order derivatives. A posteriori tests using the Taylor series approximation in LES of forced 
isotropic turbulence and channel flow confirm that the Taylor series approximation yields accu- 
rate results for the dynamic coefficient. Moreover, the simulations are stable and yield accurate 
resolved velocity statistics. 



1. Introduction 



In large-eddy simulation (LES), the Navier-Stokes equations are filtered in an attempt to isolate the 
large-scale motion in a turbulent flow. The filtered equations contain the divergence of the sub-grid 
scale (SGS) stress, 
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where (~) denotes a filtering operation. This filter, called the grid filter, is a convolution with the kernel 

fix) = / f{x')GA{x,x')dx', (2) 



which has a smoothing action on the scales smaller than A. Analogous to the closure problem related 
to the Reynolds-averaged Navier-Stokes equations, the SGS stress must be modeled using available 



filtered quantities. A variety of SGS models exist (Piomelli, 1999; Meneveau and Katz, 2000), but the 



results largely depend upon the choice of model coefficients, and the latter often have to b e tuned from 
one flow regime to another. The dynamic procedure introduced by Germano et al. (1991 ) avoids such 
ad-hoc tuning. A crucial step in the dynamic procedure is using a filtering operation, called test filtering, 
to gather information about the smallest resolved scales. The test filtering operation, denoted by (~), 
is defined in a similar way to the grid filter, except the filter acts on a larger scale aA: 



fix) 



f{x')GaA{x,x')dx', (a > 1). 



(3) 



In this paper, we restrict attention to the dynamic formulation ( Germano et al, 1991 ) of the Smagorin- 
sky ( ^magorinsky, 1963 ) eddy- viscosity model. It approximates the deviatoric part of the SGS stress 

by 

1 



-2{c,Ay\S\S, 



(4) 



where Cg is the Smagorinsky constant, Sij = ^{djUi + diUj) is the filtered rate of strain tensor, and \S\ — 
(2 Sij Sij ) is the filtered strain-rate magnitude. The expression for , obtained by minimizing the mean 



square error in the Germano identity (see Germano et al., 1991; Lilly, 1992), £ = 



L 



f 



IS 



{4f 



where the angle brackets denote an averaging operation (Ghosal et al, 1995), 



Li^ = UiUj — UiU 



is the resolved turbulent stress, and 



M . 



2A' 



\S\S,,-a^\S\S,^ 



(5) 



(6) 



(7) 



(the superscript "/" in both tensors refers to "filtering" — to be later contrasted to Taylor series ap- 
proximations). In (^, for simplicity, we have put a in place of the ratio of the 'compound' filter length 
to the grid filter length, where 'compound' filter length is the effective length scale of the filter obtained 
by sequentially applying the grid and test filters. The error in doing so is tolerable for typical values 
of a and one can calculate the precise form of this ratio after choosing a specific type of test filter and 



assuming a form for the implicit grid filter. The reader is referred to Winckclmans et al. (1998) for 
further discussion. The averaging in (0) may be done over directions of statistical homogeneity (Ger- 



mano et al., 1991), if any exist. In complex geometries without homogeneous directions, the Lagrangian 



dynamic model (Meneveau et al, 1996), which calculates time averages along pathlines, can be used. 
In the present work, only turbulence with statistically homogeneous directions will be considered for 
simplicity. Hence, spatial averaging is employed in all applications. 

However, from a practical perspective, test filtering adds computational cost to LES. The cost is 
typically manageable when dealing with pseudo-spectral numerical methods, where filter operations can 
be performed in Fourier space. When using physical-space based test filtering approaches (e.g. in finite- 
difference or finite- volume codes with structured grids), the operation co unt depends upon the number 
of neighboring grid-points involved in the filtering. Najjar and Tafti (1996 ) give a discussion of the effects 
of using test filters with finite-difference approximations and implications for the dynamic Smagorinsky 
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model. When dealing with complex-geometry flows, unstructured grids are often employed for which one 
must decide which neighboring nodes are involved in the filtering. Challenges also arise when seeking 
parallelization. These difficulties have somewhat limited the applicability of the dynamic model to LES 
of complex-geometry fiows. Various filtering operators for unstructured grids have already been proposed 
and tested, in several papers by Jansen (1994, 1999). He discusses and compares several options, including 
derivative based filtering, and generalized top-hat filtering (Jansen, 1999). The generalized top-hat filter 
is a natural extension of the top-hat filter to unstructured grids by averaging over all elements that 
share a particular node. Filtering methods for complex geometries have also been discussed by Mullen| 
and Fischer (1999|) . 

The derivative-based filter approximates the function to be test filtered by expanding it locally in 
a Taylor series and truncating. This method is attractive because typical codes already have efficient 
methods in place to evaluate de rivatives (as opposed to filte ring, which is not typically needed in most 
existing codes). As reviewed in Meneveau and Katz (2000| ), the idea of expanding the local velocity 
field in a Taylor series has been used before in SGS modeling, mainly to simplify modeling terms for 



the similarity model and/or the Leonard stresses (Leonard, 1974; Clark et al, 1979; Liu et al, 1994 



Winckelmans et i 



1 



, and also in the context of so-called defiltering SGS models ( Stolz and Adams 

1999 ; Kuerten et al, 1999| ) . In the context of the dynamic Smagorinsky model, Taylor series expansion 
was originally suggested by Gao and O'Brien (1993| ) to analytically study the limit a — > 1, but was not 
applied or extended in LES. As me ntioned above , the derivative based filter (or Taylor series approach) 
was one of the options tested by Jansen (1999), although no detailed results are presented for this 
approach. Jansen has used generalized top-hat filtering in his dynamic LES of airfoil flow because this 
method has been found to be cheaper than the Taylor series method and yet gives similar accuracy 
(Jansen, 1999), in the context of his specific code. The derivative based method has also been used to 
perform a priori studies of the sensitivity of various SGS models to the type of test filtering employed 
(Sagaut and Grohcns, 1999). 

In the present paper we focus on the Taylor series method because it is a fairly general approach 
which can be formulated, studied, and tested in general, less code-specific terms. Section ^ describes 
the formulation of the approach, in which the dynamic Smagorinsky coefficient is expressed entirely in 
terms of derivatives of the resolved velocity field and no test filtering is required. The accuracy of this 
approach is tested both a priori (§ |^) as well as a posteriori in LES of two benchmark flows, namely 
forced isotropic turbulence and minimal channel flow (§ ^). Structured grids are used both for the direct 
numerical simulation (DNS) used in the a priori tests and in the LES runs, in order to allow us to 
make comparisons in a highly controlled and standard numerical environment. Basic conclusions are 
presented in § a. 



2. Formulation 

Given a quantity f{x) (time dependence is implicit) that is to be test flltered, f{x') in Equation 
(^ (written for a homogeneous fllter) is replaced with its Taylor series expansion about the point 
X. This leads to simple integrations that are performed analytically. The result is (assuming uniform 
convergence) 

— °° f— Li" ~i 

f(^)-^T.—^^n■■■^J (2/n---2/»J, y = x-x', (8) 

^ — ' Til X 

n=0 

where now the angle brackets denote the operation of taking the mean with respect to a homogeneous 
filter GaA, e. g., 

iVii ■•■ViJ ^ J y»i • • • 2/j:„ GqzI (y ) dy . (9) 

For an isotropic filter, we write GaA{y) — GaA{\y\)^ and all terms with n odd vanish, so that we can 
take n — 2m. The choice of an isotropic filter also effects the remaining terms since they all become 
isotropic tensors. Since derivatives can only be calculated to a limited accuracy, and we are forced to 
truncate the Taylor series, this method yields an approximation to the filtering operation. By varying 
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the number of terms used in the Taylor series, and the way derivatives are calculated, the accuracy of 
this approximation can be varied. 

Here, a formulation based entirely on the Gaussian filter, 

3/2 



G 



Gauss 



(y) 



6 



exp 



(10) 



is given although the approach can be extended to any isotropic filter with finite moments (this excludes 
the spectral cutoff filter, which has infinite second moments). Additionally, up to second order only, 
the following is also valid for the box filter, since the second moments of the Gaussian and box filters 
are the same. This means that the rest of the development in this section, up to second order only, can 
also be applied when using the box filter. This fact is used later in §4.2, where the box filter is used 
to perform test filtering. Restricting attention to the Gaussian filter, the moments in (^ can easily be 
evaluated as 

where the sum is over all (2m)!/ (2'"m!) ways of partitioning {ii, i2, 
Using the result 



we have 



(aZ\)2" 



12'^ 



12 ' 



■S,. 



Using (O) in (||), we have 



m=0 



24™m! 



(V^)"/(^), 



, i2m} into pairs (Isserlis, 1918). 

(12) 

(13) 

(14) 



where (V^)"* denotes m applications of the Laplacian operator. This expression allows one to replace 
test filtered quantities that appear in the dynamic model by expressions involving resolved (grid-filtered) 
quantities and their derivatives. 

Application of (|lj) to Ui, Uj, and UiUj yields the approximation 



12 dxk dxk 



288 



dui d f d^Uj 



dxk dxk V dx"^. 



dxkdxe dxkdxe 







duj_' 


dxk 




dxk_ 



(15) 

Here, the terms involving derivatives of the product UiUj have been expanded to take advantage of the 
favorable cancellation. The superscript "t" stands for "Taylor series approximation". Similarly, after 
applying (^4|) to Sij and |S'|S'y, we can calculate an approximation for My, 

(aA)' 



Mlj 2A' 



24 



1152 



OiA"^) 



where 15*1 is the derivative-based approximation to \S\, 



15*1 = 



2 5,; 



24 



8^ 9 



1152 ^ 



0(zi6) 



1/2 



(16) 



(17) 



Here, the derivatives of products have not been expanded, since no cancellation in the terms of M*^- 
occurs. The dynamic coefficient c* can then be computed as before by evaluating tensor contractions 
and averaging over homogeneous directions. For example, below we will make frequent use of the case 
where only the 2"^^ order expansions are kept, so that we have 



{dkUidkUj) (^{\ 



dl (|, 



24 



((151 



24 



dl,(\S\S,n)-a'^\S^\dl,S,^ 



(18) 
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with 



15* 



2 Sr 



2^ '-'m'-'tj 



1/2 



(19) 



In the case of isotropic turbulence, the test filtering is done in three dimensions, so that the values of 
the indices k and m in ([l8| ) are 1, 2, 3, and the averaging in the numerator and denominator is volume 
averaging. However, in the case of channel flow, in order to approximate planar (in x-z planes) instead 
of volumetric filtering, the Laplacian term in the expansion of (^ only contains derivatives in the x 
and z directions (i. e. the indices k and m in (|l|) and (|l|) only cycle over the two values fc = 1, 3 and 
rn = 1,3). Additionally, in the channel flow case, plane averaging of the numerator and denominator of 
(|l8|) is performed in the x-z planes. 



3. A priori tests 
3.1. Methods 

The a priori tests are done using a 128^ DNS database of forced isotropic turbulence at R\ « 94, 



produced using the same pseudo-spectral algorithm as in Cerutti and Meneveau (1998). The grid filter 
size. A, is chosen to be the grid spacing of a 32'^ grid, and corresponds to A/rj ~ 9, where rj is the 
Kolmogorov scale. To perform the data analysis on the 32'^ grid, the 128^ DNS data is grid filtered (using 
the Gaussian filter) and stored back on the DNS grid. Then, the grid-filtered velocity is filtered using 
a spectral cutoff filter of width A (the same as the 32^ grid spacing) so that no aliasing occurs when 
transferring the velocity to the 32"^ grid. The grid-filtered velocity, on the coarse grid, is then used to 
calculate the filtered rate of strain tensor using finite differences. The model coefficient is calculated 
using the approximations (15), ([l^), and (17), and is compared to values obtained by explicitly applying 



a test filter. The series (|lj) must be truncated, so here two levels of accuracy are considered: second 
order and fourth order. The second order approximation (denoted as T = 2) is obtained by keeping 
terms 0{A'^) and lower in ([T^), so only the first term in (|l5|), and only the first two terms inside the 



curly braces in (16) are kept. The fourth order approximation (denoted as T = 4) is obtained in an 
analogous way and contains all of the explicitly shown terms in (|l^) and (^6|) . 

At both levels of accuracy, all of the derivatives required to implement the derivative-based test 
filtering arc initially calculated using centered finite differences accurate to second order (denoted as 
D = 2). These derivatives could be calculated spectrally, however the goal here is to examine how 
the derivative-based method performs when derivatives are calculated as they would be in a complex 
geometry situation. In an attempt to separate the errors due to truncating the series (|l|) and finite 
differencing errors, fourth order finite differences (denoted as D = A) are also tried for calculating 
derivatives. This allows one to keep the filter size constant while the accuracy of the finite differences 
is changed. 

A coarse 32'^ grid is used to store the grid-filtered velocity during the a priori test, so that all 
subsequent calculations would be the same as if a LES was actually being performed. This is also 
consistent with the choice of grid filter size. The ratio of the widths associated with the test filter and 
the grid filter, a, is varied between 2, 3, and 4, though in this study, the most attention is given to the 
common choice a = 2. 

Finally, the approximations (^5|) and ( |l6| ) are calculated to the desired level of accuracy and the 
results are used in (||) to calculate c^. Spatial averaging is employed in (||) because of the spatial 
homogeneity of the velocity field under consideration. Note that since (^5|) and ( |l6| ) contain products 
of velocity components, they also contain aliasing errors. However, these were not removed since it is 
assumed that methods to do so may not be available during a typical LES of a complex geometry flow. 
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3.2. Results 



The results obtained from using the Taylor series approximations are compared to those obtained by 
performing the classical test filtering. The comparison is based on scatter plots of tensor element values 
for some representative cases, and on correlation coefficients and normalized mean square errors among 
individual tensor elements. The correlation coefficient for a given tensor element is defined as 



\ mn mn / \ mn/ \ mn/ 



Tt 2 

mil 



1/2 



(20) 



no summation over indices), and similarly for the My tensor elements. In order to also quantify the 
agreement among the magnitudes of the tensors, we compute the normalized square error ( Liu et aT] 
l999| ) defined as 



(21) 



and similarly for the My . Also the dynamic coefficient is evaluated by averaging the tensor contractions, 
and the relative error in ci is obtained. 



3.2.1. Second order approximation 

In this section, the results obtained by using the second order approximation (T = 2) and second order 
finite differencing [D — 2) are presented. Here, the parameter a is fixed at the common value of 2; effects 



of varying a are discussed in § 3.2.5. First, we compare the results for the components Lij obtained using 
the two methods. Scatter plots comparing the 2"^* order approximation to the exact (filter-based) results 
are shown in Figure |^ for a diagonal and an off-diagonal component of Lij . Both plots indicate that the 
2nd order approximation is overall quite good, but has a tendency to under-estimate the magnitude of 
the L-components for small \Lij\, but over-estimates the magnitude at larger |Ly |. 

As shown in the first column of Table |l|, the 2"^ order approximation is well correlated (p ^ 0.95) 
with the exact results, and the normalized square error between the two is 6% for diagonal terms 
and ~ 15% for off-diagonal terms. 

Scatter plots of components of My- obtained using the 2"*^ order approximation and the exact results 
are shown in Figure |^, from which it can be seen that the approximation is very good. As shown in 
Table 1^, the 2'^'* order approximation for the My is very highly correlated with the exact values, with 
p ~ 0.99. The normalized mean square is significantly lower than it was for the Ly , being ^ 1% for 
diagonal terms and ~ 0.5% for off-diagonal terms. 

Using the approximate and exact results for Ly and Mijiii Equation (^), and using volume averaging 
gives the values of shown in the last entries of Table [l[ The relative error in obtained by using 
the 2"*^ order approximation is ~ 5%. Since the value of the coefficient obtained in the present work is 
essentially the same as for the standard filtering approach, the dissipation characteristics will be the 
same as for the standard filter-based dynamic model. 



3.2.2. Fourth order approximation 

To answer the question of whether one can improve on the 5% relative error of the 2'^'^ order approximation 
by including higher order terms, a priori tests are done with the 4*"^ order approximation (T = 4, but 
still = 2). As in the 2"*^ order case above, we fix a = 2 here. Beginning with the components of 
Lij , the scatter plots in Figure |^ show that the 4**^ order approximation tends to under-estimate the 
magnitude of the components of Lij. Despite this, the correlation coefficient between the 4"^ order and 
exact results, shown in the second column of Table |^ to be 0.98, is higher than it was between the 
2nd order and exact results. Additionally, the normalized square error between the approximation and 
exact results for the Ly , also in Table 0, is lower overall with an improvement in the off-diagonal terms 
outweighing a slight increase in the error of the diagonal terms. Note that with the use of the 4*** order 
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approximation, there is no longer a large discrepancy between the diagonal and ofF-diagonal terms, as 
was observed when the 2"^^ order approximation was used, as the normalized square error is ^ 7%. 

For Mij , Figure ^ shows scatter plots obtained with the 4*'' order approximation, which now shows 
a slight tendency to over-estimate the magnitude of Afy . The correlation coefficient is seen in table |^ 
to remain high, at ^ 0.99. However, the normalized square error for the Mij components is increased 
significantly to 3% for both diagonal and off-diagonal terms by using the 4*'^ order approximation. 
This is an increase from 1% for diagonal terms and 0.5% from the off-diagonal terms, which were 
obtained with the 2"^^ order approximation. As a result, the value of obtained with the 4*'' order 
approximation has a relative error of ~ 27%. This is significantly worse than the results obtained using 
the 2"^^ order approximation. 

There are two likely reasons why the 4**^ order approximation gives inferior results. The first is 
that the finite difference errors in the 2"*^ order terms are swamping the 4**^ order correction. In other 
words, due to the second-order finite differencing {D — 2), the error in dkUidkUj in Equation ^ is 
0{A^) — the same order as the 4*'^ order correction. The second possible reason is that the filtered 
velocity field is simply not smooth enough to allow reliable calculation of high order derivatives through 
finite differencing. This issue is an inherent, well-known difficulty associated with LES velocity fields 
that have a Kolmogorov energy spectrum. Even after removing some energy near the grid-scale through 
the implicit grid filter, gradients are still dominated by modes very near the scale A. Also note that 
the main problems seem to be associated with the error in the 4**^ order My terms, which contain the 
highest order velocity derivatives. These two possibilities are examined in more detail below. 



3.2.3. Higher order derivative scheme 

To check the possibility of whether the truncation error from the 0{A^) finite differences is swamping 
the 4"^ order correction, more a priori tests were done using 0(Z\^) finite differences to evaluate the 
2nd order terms. This ensures that the finite difference truncation error from the 2"*^ order terms is 
now of higher order than the 4'^ order correction terms. To compare the relative importance of finite 
differencing error and the error associated with truncation in the Taylor series method, we consider the 
form of the leading error terms. This is done in one dimension only for simplicity. In the case T = D ^ 2, 
the expansion with finite differences is 

7w = /w + ^(rw),,+oM^(5^^) {i"i.,)^^+o(A% (22) 

where the subscript fd refers to derivatives calculated with centered finite differences, e.g. (^f"{x)j = 

{f{x + A)-2f{x) + f{x-A))/A'^. In the case T = 4, D = 2, the leading error term is the 2/1152 part 
of the coefficient of /™ in the above expression, which is of the same order (in A) as the fourth order 
term in the approximation to the test filter (i.e. the q;^/1152 part of the coefficient of /™). In particular, 
when a = 2 the finite differencing error is smaller but comparable to the fourth order correction term 
in the Taylor series expansion for the test filtered quantities. 
With Z? = 4, the expansion becomes 

In the case T — 2, D — A, we see that the leading error term is simply the fourth order correction term 
in the Taylor series expansion of the test filtering operator (the (aZ\)^/1152 term). While in the case 
T — D — A, the leading error term is of higher order than the fourth order correction in the expansion 
of the test filtering operator (it is oc A^). When a = 2, the finite difference error is larger, though of 
the same order of magnitude as the correction term for the test filter expansion. 

The results for the and Mij components, again for a = 2, are summarized in the third and fourth 
columns of Table 0. For the 2'^'^ order approximation, although the correlation coefficients remain nearly 
the same as when the 2'^'* order finite differences were used, the normalized square error for shows 
a significant increase, while it shows a smaller increase for M^j. The resulting error in for the 2'^'* 
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order approximation with 4'^ order finite differences is mucfi larger (37%) than when 2"*^ order finite 
differences are used. When the 4**^ order correction terms are included, evaluated with 4**^ order finite 
differences, the correlation coefficients for the Lij and Mij components do not change significantly. 
In contrast, the normalized square error for the Lij obtained with the 4'^ order approximation and 
^th Qj-^Qj. differences is decreased significantly from both the 2""^ order case with the 4*^ order finite 
differences and the 4*^ order case with only 2"*^ order finite differences. The normalized square errors for 
the Mij calculated using 4*^ order finite differences are increased from both the 2"*^ order approximation 
with 4**^ order finite differences, and the 4*'^ order approximation with 2'^'^ order finite differences. This 
shows that use of 0(zi*) finite differences can reduce error in using the 4"^ order approximation, but 
overall the results are not as good as when just 2"*^ order finite differences are used to implement the 
2nd Qpder approximation. Since the leading error term in (23) is at higher than fourth order, the above 
analysis of the error terms does not clearly explain the worsening of results in going from T = D = 2 to 
T = D = A. A more intuiti ve an alysis based on the roughness of the underlying velocity fields is given 
in the following section (§ ^.2.4 ). 



3.2.4- Effects of smoothness 

As a check to see under what conditions the 4**^ order approximation does give better results than the 
2nd order approximation, the a priori tests are repeated with velocity fields of varying smoothness. 
A smoother velocity field allows more accurate calculation of high order derivatives, and also yields 
smaller errors in truncating the Taylor series. As a simple way to adjust the smoothness of the velocity 
field, an additional Gaussian filter is applied to the DNS data before performing the a priori tests. The 
width of this additional filter, or prefilter, is varied between A and 8Z\. This analysis is done here only 
to illustrate trends and is not proposed as a practical method in a simulation where we wish to avoid 
the need for filtering in the first place. 

The results obtained with the prefilter width set to A A are summarized in the fifth column of Table ^ 
When 2"*^ order finite differences are used, it is seen that the normalized mean square error is increased 
for some components by going from the 2"'^ order approximation to the 4*'' order approximation. The 
relative error in is still higher for the 4**^ order approximation. The results obtained by using both 
^th Qj-jgj. finite differences and prefiltering are shown in the sixth column of Table |l]. In this case, the 
normalized mean square error for the Lij components is decreased, while the results are mixed for the 
Mij. The final value of c^, as seen in Tabic 0in this case does show an improvement in going from the 
2nd order approximation to the 4**^ order approximation. This trend becomes more pronounced as the 
size of the prefilter is increased further. It is important to note that the resulting coefficients should not 
be compared to those obtained without prefiltering, because the velocity field in this case is actually 
different (it has been artificially smoothed out by the prefilter). 

We conclude that only when the velocity field is very smooth does the use of the 4'^ order approxi- 
mation yield more accurate results. However, this is not typically the case as LES fields are inherently 
somewhat rough. This result suggests that under typical circumstances, the 4"^ order approximation 
will yield poor results when compared to the 2"^^ order approximation. 



3.2.5. Effect of varying a 

To examine how the size of the test filter affects the accuracy of the derivative based method, a priori 
tests are done with a = 3 and 4, in addition to the base case of a = 2 discussed above. The derivatives in 
these tests are calculated using 2"^^ order finite differences and keeping the truncation to the second-order 
(i. e. D = 2 and T = 2). Since the derivative based method is based on truncating a Taylor series, one 
expects that the method will give better results for smaller values of a, for a given A. The correlation 
coefficients and normalized square errors for Ly and Mij are shown in the last two columns of Table 0. 
For a = 3, the correlation coefficients for Lij and Mij are decreased from their values obtained with 
a = 2. The normalized square error of the obtained for the 2"'^ order approximation with a = 3 are 
increased dramatically over those obtained with a = 2, especially the off-diagonal terms, which show 
errors ~ 120%. For the Mij, the error obtained with the 2"^^ order approximation also shows a large 
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increase in error compared to the case a = 2. The final results for are given in Table ^ and they show 
that the case a = 3 has a 6.41% relative error. Considering the errors in the individual components of 
Lij and M^, this result is surprisingly good, but is most likely due to fortuitous cancellation of errors. 

As seen in the last column in Table for a = 4 the correlation coefficients for the Lij are decreased 
from the a = 2 and a = 3 cases. The correlation coefficients for the My are also decreased when 
choosing a — A. The mean square error for both the Lij and the is increased significantly in the 
a — A case. The error for obtained with a = 4 is seen to be very large. Analysis of the higher-order 
case (T = 4, not shown) yields even larger errors. These a priori tests show that the Taylor series 
approximation to test filtering yields unsatisfactory results for a > 3 in this isotropic turbulent flow. 



4. A posteriori tests 

As a complement to the results of the a priori tests, a posteriori tests are used to compare the Taylor 
series based dynamic model to the traditional method based on test filtering. The tests are done by 
performing LES of both forced isotropic turbulence and channel flow. Since the best results in the a 
priori tests were obtained by truncating the Taylor series expansions at second order, all a posteriori 
tests are done using T = 2 and D = 2. 



4.1. LES of Forced Isotropic Turbulence 

To test the derivative-based method in isotropic turbulence, the pseudo-spectral DNS code is modified 
to perform LES on a 64^ grid. This includes adding a 3/2-rule deahasing procedure, which is apphed 
only to the nonlinear term in the filtered Navier-Stokes equations {u x lD) and not the SGS stress 
term. This is done because the nonlinear term is known to be exact, while the SGS term is already an 
approximation. In addition, the type of nonlinearity contained in the SGS stress term is not completely 
removed by standard dealiasing procedures such as t he 3/2-rulc or random phase shifts. Simulations 
are forced with a constant energy injection rate (as in Cerutti and Meneveau, 1998| ) of e = 0.0007, the 



molecular viscosity is v = 10~ , and the domain is of length L = 27r. At statistical steady state, the 
Taylor-scale Reynolds number of the flow is R\ « 2100. Simulations using both exact test filtering and 
the Taylor series based method are performed. As with the a priori tests, a Gaussian test filter is used 
at scale 2 A (corresponding to a = 2). The exact test filtering is done in Fourier space, while the Taylor 
series method uses second order finite differences to calculate all of the derivatives associated with the 
dynamic procedure. The initial condition is a random field with a k~^^^ spectrum and random phases. 
Volume averaging is used in this spatially homogeneous flow. 

Figure ||(a) shows the time evolution of Cs for exact test filtering and the derivative based approach. 
Although the derivative based method consistently gives a value of Cs roughly 5 % higher than the 
exact method (in contrast to the a-priori tests), both methods relax from the initial value (cs(0) = 0, 
as expected for Gaussian fields) to their quasi-steady value at about the same rate. Most importantly, 
the calculation done using the derivative based method is seen to be stable. The quasi-steady values for 
both methods lie near Cs ~ 0.13 and are lower than the value of Cs ~ 0.16 often observed in isotropic 
turbulence. This may be a result of using a Gaussian filter. As shown in Figure |^(b), the radial energy 
spectra obtained from the two methods are essentially the same, only showing a slight difference at 
high wavenumber due to the slightly different values of c^. Both cases agree well with the Kolmogorov 
—5/3 prediction, with ck = 1.6. 



4.2. Application to channel flow 



To examine the effects of anisotropy and the presence of walls on the derivative-based procedure, 
moderate Reynol ds number cha nnel flow simulations are performed with a LES version of the DNS 
code NTMIX3D ( Stoessel, 1995 ). In this code for fully compressible flow, the governing equations are 
integrated using an explicit third order Runge Kutta time advancement and a sixth order compact 
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finite difference scfieme (Leie, 1992). No siip isotliermai boundary conditions are used at the top and 
bottom boundaries, whiie periodic boundary conditions are used in the streamwise (x) and spanwise 
(z) directions (which are homogeneous directions). A driving source term compensates the shear stress 
at the walis which ailows to run the simuiations using periodicity in the x direction. A computational 
domain of size — 27r, Ly = 2 {= 25), and Lz = 0.908 in the three directions is used. The x, y, and 
z coordinates belong to the intervals < a: < L^, —Ly/2 < y < Ly/2, and < z < L^- 

The shear Reynolds number (based on friction velocity Ur and channel half width S) is Rcr = 180, 
which corresponds to a convective Reynolds number Rbc — 3300 (based on channel half width and axial 
velocity). To allow a reasonable time step within the restrictions of the acoustic CFL condition, the 
mean centerline Mach number M — 0.2 has been chosen in the low subsonic domain. The trace of the 
sub-grid scale stress tensor can be rewritten as Tkk = I^sgsP- small Mach number considered 

in the present simulations, the sub-grid Mach number Msgs is expected to be small. Consequently 
we simply neglect Tkk- Also, SGS fluxes in the total energy equation were treated using a fixed (non- 
dynamic) value Prf = 0.6. Because of the low Mach number, the linkage to the momentum field was 
totally negligible. The simulation parameters have been fixed to allow quantitative comparisons with the 
direct numerical simulation data of Kim et al. (1987 ). Two LES runs with the dynamical Smagorinsky 
model are performed: one with the classical filtering procedure and the other with the Taylor series, 
derivative-based method. 

For both LES cases the computational grid contains 17 x 61 x 16 points in the x, y, z directions 
(or xi, X2, X3). The grid is uniform in the streamwise and spanwise directions and the corresponding 
resolution was A+ w 66 and A+ w 10. Following |Gamet (199"9| ), we use a non-uniform mesh in 
the wall normal direction based on the distribution: yi — -^ta.iih{K r]i) with K = atanh(^) and 
— 1 < rji — 2^ - 1 < +1. The constant C is such that zi+ « 2 at the wall and Z\+ « 10 near 
the centerline. For comparison, minimal channel DNS was performed using = 34, Ny — 121 and 
Nz — 32. LES are started from filtering a fully turbulent DNS field. The statistics were accumulated 
over time tUr/S = 5.4. 

As usual in channel flow simulation with the dynamic procedure, we apply the test filter operation 

only in the streamwise and spanwise directions. Because of the non uniform mesh in the cross direction, 

1/2 

we choose to use for the test filter width: A — a (A^Az) ' , with a ~ 2. For the LES of the filter-based 
reference case we apply a box filter, using a trapezoidal rule for the integration of the convolution 
operation. For the derivative-based method, the sec ond order approximation [T — 2, D — 2) has been 
chosen based on the results presented in section 3^. Here the equivalence at second order between the 
Gaussian and box filters mentioned in §2 is used. The final expression for (c*) is given in (|l8|), with 
planar Laplacians {k = 1,3 and m = 1,3) and planar averaging, as discussed in § 0. 



4.2.1. A priori results based on initial condition for LES 

As a first step we perform a priori tests at the initial time of the LES simulations, which are filtered DNS 
fields. Scatter plots comparing the Ln, L12, Mn, M12 components obtained by the Taylor series and 
filtering based methods, in the x-z plane at j/"*" = 43, are shown in Figure |^. The agreement between 
real filtered and estimated local values is fair. The general trends are reproduced with correlation 
coefficients of order 0.7. Specifically p(Ll^, = 0.79, piL\2,L{^) = 0.81, p{M[^,Ml^) = 0.74, and 
p{M[2, ^12) = 0.64. Even if the local behavior of the typical components of Lij and Mij is not highly 
accurate, the global behavior of these components is estimated reasonably well, as demonstrated by 
Figure 0, showing the same tensor elements averaged in the x-z planes, as a function of wall normal 
direction. We notice that despite an overestimation of the magnitude of each component, their general 
shape is correctly predicted by the derivative-based approximation. 

Figure Ijshows the instantaneous Smagorinsky coefficients obtained by the classical filtering procedure 
and by the present derivative-based method. The a priori test shows that the Smagorinsky coefficient 
obtained by the Taylor series expansion approach is slightly smaller than the classical filter-based results. 
However, the agreement is sufficiently good in the context of a practical scheme, especially the reduction 
to zero in the well-resolved near-wall region. As we will see in the next section, a posteriori tests give 
improved results. 
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4-. 2. 2. A posteriori results oj LES 

All the results presented in this section are time-averaged over 55 realizations spaced by a time interval 
of 0.16 /ur- In Figure ^(a), we observe that the Smagorinsky coefficient is very well predicted by the 
derivative-based method, except some deviation near the centerline. As shown in Figure |(b) the 
derivative-based approach yields the expected y"*"^ behavior in the vicinity of the wall. 

In Figure ^0|, we observe a fairly good agreement in streamwise mean velocity U between the results 
of the Taylor series expansion method and both the results of the classical filter-based dynamic model 
and the DNS results. Despite a slight underestimation in the buffer region, the proposed approach leads 
to quite accurate results. 

In Figurejl^ are shown the rms intensities of the large-scale (resolved) turbulent velocity fluctuations. 
These rms quantities are compared with the rms velocities from the DNS filtered at scale A using a box 
filter. The agreement between both LES results and the filtered DNS data is quite good. For the derivative 
based method, the streamwise rms velocity is slightly overestimated by both LES approaches. The wall- 
normal and transverse rms velocities are underestimated slightly more by the derivative-based method 
than the filter-based method. In Figure |ll|(d), we compare the Reynolds shear stress of the resolved 
velocity, {u'v'). The proposed approach yields very good agreement with the filter-based dynamic model 
and the reference filtered DNS data. 



5. Conclusions 

An implementation of the dynamic Smagorinsky model that uses Taylor series expansions to avoid 
explicit test-filter operations in LES has been developed. The proposed method has been subjected to a 
priori and a posteriori tests. These tests were performed using structured grids, to check the performance 
of the method in idealized, well controlled reference cases. 

The a priori results obtained by truncating the Taylor series at 2"^^ order and at 4"^ order have been 
compared with results obtained by test filtering. It was found that for LES of isotropic turbulence, it is 
possible to obtain values of accurate to 5 % by using the 2"^^ order approximation. However, results 
are not improved by using the 4*^ order approximation because of the errors associated with evaluating 
high-order derivatives on inherently rough LES fields. From these a priori tests, it is concluded that 
the derivative based method should be implemented using the 2"^^ order approximations, with a ~ 2. 
A posteriori tests of forced isotropic turbulence show that the derivative-based method yields stable 
results, and values of to within 10% of those obtained by explicit test ffltering. 

Applications to a minimal channel configuration show that the observed differences between the 
tensor elements in the filter-based and Taylor series based approaches are larger than those in isotropic 
turbulence. This could be due to the strong anisotropy of the test filter as well as of the turbulence 
itself. However the Smagorinsky coefficient is correctly estimated by the proposed approach, especially 
the y^^ behavior in the vicinity of the wall. Moreover ffist and second order statistics, which are the 
most important ones for practical engineering calculations, are correctly predicted when using the 
derivative-based dynamic model. 

Strictly speaking, these conclusions are applicable only to the fairly simple flow configurations 
considered in this work. However, the results suggest that applications of the Taylor series based dynamic 
model to co mplex-geometr y flows on unstructured grids is a promising direction. Taken together with 



the results of Jansen (1999 ) on the alternative generalized box- filter approach, the present results show 



that applications of the dynamic model to LES of complex-geometry turbulent flows are feasible. 
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Table 1. Correlation coefficients, normalized square error, values of c^, and relative error in obtained from a priori tests 
in forced isotropic turbulence. 



prefiltering 



order of finite dif- 
ferences {D — ) 



order of 
sion (T 



expan- 



p(i(i,Mi) 

p{Ml^,Mi,) 
p{M{^,MU) 

^(-^'121 -^12) 

f(M/i,Mi*i) 
£{M(^,Mi2) 

{cif 

error in (cl)^ 



(%) 



4A 



0.955 
0.947 
0.996 
0.998 

0.056 
0.158 
0.009 
0.006 

0.0151 
0.0159 
5.30 



0.978 
0.976 
0.990 
0.993 

0.070 

0.078 
0.031 
0.025 

0.0151 
0.0110 
27.2 



0.946 
0.932 
0.995 
0.996 

0.187 
0.514 
0.014 
0.009 

0.0126 
0.0173 
37.3 



0.981 
0.977 
0.986 
0.989 

0.029 
0.050 
0.047 
0.040 

0.0126 
0.0104 
17.5 



0.993 
0.993 
1.000 
1.000 

0.023 
0.048 
1.08E-4 
1.02E-4 

0.0097 
0.0110 
13.4 



1.000 
1.000 
1.000 
1.000 

5.78E-4 
9.13E-4 
3.03E-4 
4.51E-4 

0.0097 
0.0094 
3.1 



0.889 
0.863 
0.832 
0.868 

0.369 
1.264 
0.363 
0.268 

0.0156 
0.0146 
6.41 



0.791 
0.748 
0.380 
0.415 

1.301 

4.556 
10.315 
7.845 

0.0164 
0.0005 
97.0 
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(a) (b) 

Figure 1. Scatter plot of typical components of Lij evaluated using the Taylor series approach compared to the classical 
filter- based approach, for the case a = 2 with the 2'"* order approximation (T = 2). (a) shows Lii while (b) shows -Li2- 





Dynamic Model for LES Without Test Filtering: Quantifying the Accuracy of Taylor Series Approximations 



15 





(a) (b) 
Figure 4. Same as in Figure ji] using the the 4"^ order approximation (T = 4). (a) shows Mn while (b) shows M12. 
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Figure 5. Results from a posteriori tests in LES of forced isotropic turbulence on a 64=* domain with _R;^ 2100. (a) 
time evolution of dynamic Smagorinsky coefficients. Solid line: derivative-based, dashed line: test-filter based, (b) radial 
energy spectra. Solid line: derivative-based, dashed line: test-filter based dynamic model, long dashed line: normalized 
Kolmogorov spectrum cx{kr})~^i'^ with cx = 1-6. 
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Figure 10. Mean Velocity U in wall coordinates: — DNS, x derivative-based method, • classical filter-based dynamic model. 
Dashed line: U/ut = 2.b\a{y+) + 5.5. 



